function [G, V] = speciter2d(N);
% Solve the Green's function of -Delta u + (1+sin(2*pi*x)*cos(4*pi*y)) u  on 
% [0,1] * [0,1] with periodic boundary condition with spectral method
% and CG iterative scheme.
h = 1 / N;
mesh1d = (0:N-1)'*h;
[X,Y] = meshgrid(mesh1d, mesh1d);
k1d = fftshift(((-N/2):(N/2-1))');
[KX, KY] = meshgrid(k1d);
ksqrmat = 4*pi^2*(KX.^2 + KY.^2);
V = 1+sin(2*pi*X).*cos(4*pi*Y);
ksqrrep = repmat(ksqrmat, [1 1 N N]);
Vrep = repmat(V,[1 1 N N]);


G = zeros(N*N, N*N);
f = eye(N*N,N*N);
f = reshape(f,N,N,N,N);
nrhs = N*N;
maxiter = 500;
u0 = f;
r = f - specA2d(N,u0,Vrep,ksqrrep);
w = -r;
z = specA2d(N,w,Vrep,ksqrrep);
a = sum(sum(r.*w,1),2) ./ sum(sum(w.*z,1),2);
a = repmat(a,[N N 1 1]);
ux = u0 + w .* a;
B = 0;

for iter = 1 : maxiter
  if( mod(iter,10) == 0 )
    fprintf('%5d-th iteration\n\n',iter);
  end
  r = r - z .* a;
  if( sqrt(sum(sum(sum(sum(r.^2,1),2),3),4)) / nrhs < 1e-8 )
    break;
  end
  B = sum(sum(r.*z,1),2)./sum(sum(w.*z,1),2);
  B = repmat(B,[N N 1 1]);
  w = -r + w .* B;
  z = specA2d(N,w,Vrep,ksqrrep);
  a = sum(sum(r.*w,1),2) ./ sum(sum(w.*z,1),2);
  a = repmat(a,[N N 1 1]);
  ux = ux + w .* a;
end
G = reshape(ux, N*N, N*N);
